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We investigate the nonlinear evolution of the Bales-Zangwill instability, responsible for the mean- 
dering of atomic steps on a growing vicinal surface. We develop an asymptotic method to derive, in 
the continuous limit, an evolution equation for the two-dimensional step flow. The dynamics of the 
crystal surface is greatly influenced by the anisotropy inherent to its geometry, and is characterized 
by the coarsening of undulations along the step direction and by the elastic relaxation in the mean 
slope direction. We demonstrate, using similarity arguments, that the coalescence of meanders and 
the step flow follow simple scaling laws, and deduce the exponents of the characteristic length scales 
and height amplitude. The relevance of these results to experiments is discussed. 
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I. INTRODUCTION 

The homoepitaxial growth of semiconductor and 
metallic vicinal surfaces, as revealed in experiments 
of molecular beam epitaxy, is ch aracter ized by a va- 
riety of morphological instabilitie^i'^IIIll The surface 
evolves by step flow driven by the external flux and 
controlled by surface diffusion of adsorbed atoms, at- 
tachment rates (Schwoebel barriers) and elastic interac- 
tions. The interplay of different physical mechanisms 
can lead to the formation of characteristic structures 
at the nanometer scale, in the form for example of 
macroscopic bunches of steps, or periodic meanders of 
monatomic steps. These spontaneously nanopatterned 
surfaces can be useful as templates in an experimental 
two steps-process, to obtain by subseq uent heteroepitaxy, 
self-organized arrays of quantum dot d^ l ^ l ^ l The physical 
description of the epitaxial surface growth depends on 
the characteristic lengths involved in the surface dynam- 
ics, from atomic processes (surface reconstruction and 
faceting) to gross macroscopic features (self-organization 
of nanost ruct ures )(212J At the mesoscopic scale of steps 
and terra ces, th e standard approach of Burton, Cabrera, 
and FranlPSIll] completed with appropriated expressions 
for the equilibrium concentration and kinetic conditions, 
allows a fairly complete description of the vicinal sur- 
face. However, the macroscopic behavior of the vicinal 
surface, in particular the longtime nonlinear evolution of 
step bunches and meanders, can be more suitably ac- 
counted by a continuum model in terms, for insta nce, o f 
partial differential equations for the surface heightl^^E^. 

In the continuum approach, initiated by Herrin^l^and 
Mullin^I^, the surface morphology is determined by the 
surface energy (75, which is anisotropic in general) and 
the surface diffusion (characterized by the diffusion coef- 
ficient Ds). It is easy to incorporate into this model other 
thermodynamical contributions as, for instance, the elas- 
tic energyS^ To be more specific, let us consider the 
epitaxial growth of a crystal under a flux F of atoms, 
relevant to molecular beam epitaxy experiments J^^ I ^^ ^^ 
We denote a the height of an atomic layer. The front 
shape is given by the function h{x,y,t) of the Cartesian 



coordinates {x, y) and the time t. The thermodynamical 
state of the system is specified by a free energy func- 
tional J-[h] of the crystal profile; the chemical potential 
per monolayer will be fj,{h) = a6J-/Sh — as in the Cahn- 
Hilliard model.E^ Therefore the conservative dynamics of 
the interface satisfies the Mullins equatior^ 



1/, = ^V72 

dt knT 



(1) 



where fee is the Boltzmann constant and T the tem- 
perature. Equation ([ij is a mass conservation equation 
with the current j ^ —W^{h). In the simplest case the 
chemical potential is proportional to the surface curva- 
ture /i = a^^sK, and using a linear approximation, one 
obtains the evolution: 



d_ 
Ft 



h = -Bs'^'^h + a^F , 



(2) 



where we defined a mobility Bs = a*Dsjs/kBT. How- 
ever, this approach, based on a free energy functional, 
valid for near-equilibrium conditions, prove to be insuffi- 
cient when kinetic processes become important. Indeed, 
the attachment of adatoms to steps modifies the mobil- 
ity coefficients (like Bg), and, coupled to external fluxes, 
can introduce new effects not related to a chemical po- 
tential. In particular, the appearance of step flow insta- 
bilities, bunching or meandering, cannot be described, in 
the continuum limit, by an evolution equation such as 
([iJ with n derived from a variational functional. One 
example is the Bales-Zangwill instabilitj^i^, which in the 
weak no nlinea r regime is described by the dimensionless 
equatioiiSSM]^ 
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dt' 



dy"^ 
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(3) 



where u — u{y, t) is the rescaled step shape {y is the 
step- wise direction). The flrst term on the right hand 
side, which is responsible for the instability, vanishes if 
the flux or the kinetic attachment barriers are absent. 
The third, nonlinear, term is also proportional to the 
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flux and cannot be derived from a free energy functional. 
We investigate in this paper, the influence of the vicinal 
surface anisotropy on the two-dimensional dynamics of 
the meandering instability. 

The meandering instability, first investigated by Bales 
and ZangwilP^, results from the difference in the attach- 
ment kinetics of adatoms approaching a step in opposite 
directions, from the upper and lower terraces. The dis- 
tinct neighborhood of an adatom sitting just above or at a 
step implies a difference in the height of the energy barri- 
ers that introduces a difference between the lower v+ and 
upper V- attachment coefficients, the so-called Ehrlich- 
Schwoebel effect)^^^. Neglecting evaporation effects the 
instability growth rate, as mentioned above, is propor- 
tional to the flux and to the difference in kinetic coeffi- 
cients a{q) ~ F{vj^ — i'-)q^, where q is the wavenumber 
of the perturbation. The original Bales-Zangwill linear 
analysis was extended to take into acco unt v arious ef- 
fects, such as desorptiorpSl kink barrierd^HH] diffusion 
anisotropyfsS or elastic interaction Although exper- 
imental measures of the instability growth rate remain to 
be performed, ther e are de tailed observations of step me- 
andering in metal^SSISSIIII^ and in semic onductors w here 
it is often associated with step bunching,'22llll2lllllll]that 
reveal the longtime nonlinear stage of the instability. 

The longtime evolution of the Bales-Zangwill insta- 
bility was thoroughly investigated in the strong nonlin- 
ear regime, by means of a continuous equation for the 
meander shape, derived from a multiscale analysis of 
the adatom diffusion model, that incorporates the effect 
of line step diffusion (related to kink barriersJP^ISHIlIIini 
and the elastic interactions between stepJ^Il ([n a one- 
dimensional approximation). This strong nonlinear 
model reproduces at least qualitatively, the meandering 
observed in metald^il when desorption and nucleation can 
be neglected. The two-dimensional effects, together with 

mound for mation were studied using kinetic Monte-Carlo 
simulationPEHmilll], 

In the case of semiconductors the coexistence of dif- 
ferent kinetic instabilities leads to two-dimensional mor- 
phologies characterized by step bunches and meanders. 
In some circumstances, bunching and meandering is also 
present in metal^SMZI Instability mechanisms operating 
in Si(OOl), based on the diffusion anisotropy that results 
from the surface reconstruction of alternating terraces, 
were found for bunchin^^ and meanderinj2fil One of 
the salient features of these two-dimensional patterns is 
their coarsening: roughening amplitude and length scale 
follow growing power laws. A unified continuum model 
of the bunching and meandering instability is lacking, 
even if phenomenological mode ls can q ualitatively ac- 
count for some of its propertie^^^EMSl The difficulty 
is to systematically derive, from the mesoscopic adatom 
diffusion model, a continuum limit that in addition to 
the kinetic processes, takes into account the step elas- 
tic interactions in its full two-dimensional form. In this 
paper we obtain such a model in the simplest physical sit- 
uation, where only the meandering instability is present 



and local step interactions are considered. Under these 
assumptions we can develop a weak amplitude nonlinear 
expansion of the Burton, Cabrera and Prank equations. 
Our two-dimensional model of the meandering instabil- 
ity contrasts to previous ones which treated the nonlinear 
evolution of an in-phase pattern, thus reducing the prob- 
lem to the (one-dimensional) dynamics of a single step. 

In the following section ^Iljwe present the basic equa- 
tions, essentially the Burton, Cabrera and Frank model. 
Section pll| deals with the weak nonlinear expansion of 
the adatom diffusion equations, that allow us to derive in 
the continuum limit a differential equation for the surface 
height. In §IV| we study the evolution of the system by 
the numerical integration of the vicinal surface equation, 
and we derive using a self-similar solution, the asymptotic 
scaling laws for the amplitude and characteristic lengths. 
The last section ^V] presents a discussion about a possi- 
ble generalization of the model to include step bunching, 
and a concluding summary. 

II. STEP FLOW: ADATOM DIFFUSION 
EQUATIONS 

The geometry of the vicinal surface schematically rep- 
resented in Fig. [l] can be described by the set of curves 
X — Xn{y,t) representing the steps n — 0,1,..., TV at 
height z = Zn = {N — n)a. The terrace T„ of initial size 
Iq, is bounded by the upper step and the lower step 
x„; its slope is then —m = —a/lo. We define the external 
normal and curvature of step n by 

_ (1, —dyXn) dyyXn , 

~ [1 + (a,X„)2]l/2 ' - [1 + (9,X„)2]3/2 ' 

respectively. 

In the absence of an external fiux F, there is an 
equilibrium concentration Ceq,n of adatoms on the sur- 
face. This concentration, which is in general not uni- 
form, depends on the curvature and elastic interactions 
of steps; it is determined by the surface chemical poten- 
tial ^iniy, t) = fln^ + fln^: 

C^e9,„ = Coe^"/^-^, (5) 

at temperature T. Cq is the adatom concentration corre- 
sponding to the reference vicinal surface, the one consist- 
ing of equidistant straight steps. The first contribution 

(c) 

/i„ ' to the chemical potential takes into account the step 
curvature, and is given by the Gibbs-Thomson relation, 

A^n ' = ^IsKn , (6) 

with 7s the step stiffness (fi = the atomic area). 
The second contribution /in^ results from the dipole mo- 
ments created by the broken bonds at each step, and de- 
pends on the distance x between steps (see for instance 
MarchenkcPl or Duport et al.l^: 

' = ^Ps{gn+1 - 9n) , ffn = 7 , (7) 
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FIG. 1: Schematic geometry of the vicinal surface, formed 
by a set of terraces r„ separated by atomic steps Sn at a; = 



where (3s — 4(1 — v'^)m1/'KElQ has the dimensions of an 
(elastic) energy per unit surface, v and E are the Poisson 
ratio and the Young modulus respectively, and is the 
force dipole moment. We assume that step n only in- 
teracts with its nearest neighbors; this approximation is 
justified in the case of the rapidly decreasing dipole 
interaction. In fact, summing over all steps would not 
change the form of the elastic repulsion in the continuum 
limit, but only renormalize the dimensional coupling con- 
stant /3s; as demonstrated by Xian^^, the ratio between 
the infinite series and one term of the sum ([7| is 7r^/6. 
Under usual experimental conditions the energies asso- 
ciated to the step stiffness and the elastic repulsion are 
much smaller than the thermal energy, so the equilibrium 
concentration at step n can be written as. 



Co [1 + ikHn + P{gn+1 - gn)] 



(8) 



where we introduced the nondimensional parameters 

7 = n-fJkBTlo , P = nps/kBT . (9) 

It is worth noting that the "equilibrium" concentration 
changes with the geometry of the vicinal surface through 
the curvature of steps and their separation, it is used in 
fact as a reference to compute the supersaturation. 

In the presence of an external flux, the evolution of 
the adatom concentration C„ = Cn{x,y,t) is well de- 
scribed by a quasi-stationary diffusion equation, as orig- 
inally proposed in the model of Burton, Cabrera and 
FrankP, 



DV^Cnix,y,t)+F = 0, 



(10) 



where D is the adatom diffusion coefficient, and the 
source term F is the deposition flux (number of atoms per 
unit time and unit surface). In ( [To| we neglected the time 
derivative by assuming a weak flux regime F <C D/^IIq. 
Under molecular beam epitaxy experimental conditions, 
with F a few monolayers per minute, this 'adiabatic' ap- 
proximation of the diffusion equation is usually satisfied. 
Also the processes of evaporation and, eventually, nucle- 
ation of surface atoms are disregarded, their characteris- 
tic times being much more longer than the typical flow 
step time l/Fl^. 



At steps, the boundaries of a terrace, adatoms are at- 
tached with a velocity rate if they come from the up- 
per terrace, or ly^ if they come from the lower one. The 
Bales-Zangwill meandering instabilitjE^ appears i n the 
case where i^^ > the Ehrlich-Schwoebel 

effectP^D, 

These parameters characterize the attachment kinetics 
that controls the flux of adatoms, then fixing their con- 
centration at the terrace boundary: 

DUn^l ■ VC„ = iy+{Cn - Ceq^n^l), X = Xn-1 , (H) 

for the upper step and, 

Drin ■ \7Cn = -i^-(C„ - Ceq.n), X = Xn , (12) 

for the lower one. The step flow results from the balance 



between the diffusion fluxes (111 and (12 1 at each step 



The normal velocity is given by 

Vn = nDnn- (VC„+i - VC„), X = Xr. 



(13) 



Equations (lOl, (111, il2l and (131 form a complete sys 



tem from which we can determine the relevant physical 
parameters. It is convenient to introduce a reduced con- 
centration, 



C-Cq 
Co 



(14) 



and nondimensional parameters related to the flux and 
the attachment coefficients, 



III 
DCi 



^±'o 



(15) 



as well as to normalize lengths with Iq (the terrace width) 
and time with Iq/UDCq, as derived from the normal ve- 
locity (13 1. Using these parameters we can write the set 



of equations in nondimensional form: 

W^Cn{x,y,t) + f = 0, 



(16) 



n„_l • Vc„ = a+[Cn - JUn-l - (3{gn - gn-l)] , 

at a: = Xn-i , (17) 

nn ■ Vc„ = -a-[c„ - 7K„ - /3(5„+i - gn)] , 

aX X — Xn , (18) 



(Vc„+i - Vc„) , a,t X = Xn- 



(19) 



Equations ( T6p9 1 describe at a mesoscopic level the dif- 
fusion of adatoms on a terrace n of a vicinal surface, 
driven by the external flux / and controlled by the at- 
tachment kinetics a±, the stiffness 7, and the elasticity 
f3 of the bounding steps. 
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III. VICINAL SURFACE EQUATIONS IN THE 
CONTINUUM LIMIT 



A stationary train of straight steps, advancing at con- 
stant velocity Xn = n + ft is an unstable solution of the 
step flow. A perturbation ^(y, n, t) of the step profile, 



Xn = n + ft + ^{y, n, t) , 



(20) 



induces a modification of the terrace adatom concentra- 
tion in the form 



Cn{x, y, t) = c„(x + ft) + (j){x, y, n, t) , 



(21) 



where the first term is the stationary parabolic concen- 
tration: 



Co (a;) 



2 



/(2 + Q_)(a+x- + l) 



solution of (16 1 with the boundaries (17 1 and (181 in the 



-1 < X < 0, 
(22) 



case of straight steps, ^(y, n, t) = 0. Note that as the flux 
/ tends to zero, the step velocity and the concentration 
vanish. 

We are interested in obtaining the growth rate of the 
meandering instability in the long wavelength, continuum 
limit. In this case it is sufficient to consider the in-phase 
mode with 



and 



{x,y,n,t) = <j>{x,n,t)e"'y , 



^(y,n,t)=e(Oe'('="+«^) 



(23) 



(24) 



where (j){x,n,t) and ^{t) are small amplitudes, and (fc,g) 
is the wavevector of the perturbation. The long wave- 
length limit is obtained by making a development in 
powers of the wavevector (up to the forth order). To 
compute the continuum limit we introduce the notation 
An = {n + 1) — n for the separation between two steps. 
In this limit both lengths, the steps separation Iq and the 
step height a tend to zero, but the slope m — q/Iq must 
be kept constant. Therefore, we formally have. 



Xn±i = n± An + ^{y,nzL An, t) , 



and assume that ^ is a smooth function of n, in the limit 
An 0. The elastic term at a: = x„ in the boundary 
conditions becomes 



and a similar expression for the boundary at a; = Xn-i, 
where the numerator in An^ ensures the correct limit 
when An (we omitted the dependence in y and t). 
We know that the meandering instability is driven by 
the flux / and the difference in attachment coefficient 
5 = q;+ — Q!_. We replace (20 1 and (21 1 in the diffusion 



equations and retain the lowest order in An, /, and S to 
get the real part of the dispersion relation: 



^q' - iq' - + mk'q' , (27) 

4a 2 2 



and the imaginary part 



2a 



(28) 



where a — a- and the amplitudes grow as ip, ^ ^ gcrt-iwt_ 
The first term in ( p7| is the destabilizing one, it is propor- 
tional to the product fS. K the long wavelength limit is 
introduced by replacing k, q ^ ek, eq, in order to keep all 
terms we must assume that, near the instability thresh- 
old, the parameters can be considered small such that: 



f^ef,S^ 7 ^ e^, (3 



(29) 



where the last two relations guarantee the balance of the 
instability growth term with the stiffness and elastic re- 
laxation terms. Using these relations all terms in a be- 
come order 0{e^), and order ©(e**) in w, suggesting the 
introduction of two time scales, e^t (instability and re- 
laxation) and e^t (dispersive waves). 

To summarize, the weak amplitude expansion of the 
step shape must be of the form. 



Xn{y, t) = n + e^{ey, n; e^t, eH) . 



(30) 



(25) and that of the concentration. 



Cn{x, y, t) = ^ e™c,„(a;, ey, n; e^t, eH) , (31) 



[An + ^(n + An) - S,{n)f 



An^ 



[An + ^(n) -^(n- An)]3 



where the variable x is itself a function of e through the 
boundary conditions, and the scaling (29 1 is assumed. 



(26) The equations of the model, become in terms of the slow 
variables: 
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dj;3^c{x, y, n) + e^dyyc{x, y,n) + ef = , 



(32) 



with 

{l,-e'^dy£,{y,n- An) ■ {dxc{x,y,n),edyc{x,y,n)) 

at X = n — An + e^{y, n — An), and 

(1, -e^dy£,{y, n) ■ {d^c{x, y, n),edyc{x, y, n)) ^ 
(l + e4[a^^(y,n)]2)i/2 



{a + e^S) 



c{x, y, n) — e^(3G^{n ~ A7i) + 



5 dyyi{y, n - An) 

^ ^ (\^i^{dy^{y,n~An)Yfl^_ 



(l + e4[a,e(y,n)]2)3/4 



at a; = n + e^(y, n), where 
a(n) = 
and the normal velocity, 



An' 



An' 



{An + e^{y, n + An) — e^(y, n))^ {An + e^(y, n) — e^(y, n — An))-^ ' 



{l,-e'dy^{y,n) 



-7^ • {da:, edy) [c{x, y,n + An) - c{x, y, n)] 



(33) 



(34) 



(35) 



(36) 



where the implicit dependence on the time variable is 
understood (time derivatives appear only in the explicit 
expression of the normal velocity). Using the expansion 
(30 31 1 together with the Taylor series in An of ^(y, n ± 
An), one may solve (|32| with the boundary conditions 



( 33p4 l to obtain a series in e of the normal velocity (36 1 



a few terms of this series are: 

Vn ^e\f (^Andr.i+An'^dr.nni 

3 

- An' -aPdnnnnS, 

+ e^An'^dnyyi 



e An \^f3 + An-a7 ) dnnyyi 

e'^An-fdyyyy^ 



e'An' 



4^9yy^+i,9yy{dyO' 



where the higher order terms (h.o.t), terms having higher 
order derivatives or powers in the step shape amplitude 
^(y, n), are shown to be irrelevant in the continuum limit. 



Note that in (37 1 the derivatives of ^{y,n) with respect 



to n are themselves functions of the small parameter e. 
This can be made explicit using the representation of the 
surface in terms of the height function h{x, y, t) satisfying 
the compatibility condition, 



J^{n; y, t) = Zn + h{xn{y, t),y,t) = 0, 



which signifies that h = const, at a height level z„ = 
{N — n)a/lo, X = Xn{y, t) = n + e^{y, n, t). In this repre- 
sentation the steps are the level contours of the surface 
z = h. 

Using the derivatives of T{n, y, t) with respect to n and 
y one obtains the derivatives of ^ in terms of derivatives 
of h. Typical derivatives are: 



dxh 



dyXn — 



9yh „ 
Oxh 



m'dxxh 



(39) 



etc. Moreover, the height itself deviates slightly from the 
mean slope m (in the continuum longwave approxima- 
h.o.t , (37) tion) h{x,y,t) — —mx + tu{tx,y,t). Keeping terms up 
to order (even number of derivatives) and (odd num- 
ber of derivatives) that correspond to the two time scales, 
one obtains the normal velocity expressed in terms of the 
reduced height u{x,y,t), with {x,y,t) the slow variables 
and with An = e = 1, 



dtu — ~dxxxxU AdxxxU 
(38) - dyy [u + dyyu + {Oyu)' + B 8 xU + C 8 xxu] (40) 



6 




y (c) y 

FIG. 2: Surface height evolution, (a-d) t= 1000,2000,3000,5350. The system size is 1024^. 



(d) 



where 



Lx, Ly (in both directions), time T and ampUtude U, 



A- 



B 



C 



/ 
S 

f_ 
5 



1/2 ^ 


2^7 ^ 






1/2 ^ 


2 






y/2 








U/3 



1/4 



1/4 



1/2 



_ (24a^/3)i/4 



2 



1/2 



T=167 



4ma7 



and where we used rescaled quantities in units of length 



and a moving frame with velocity ft (this eliminates the 
first order derivative in x). Equation ( [40| is the base of 
our numerical study of the morphological evolution of the 
vicinal surface under the Bales-Zangwill instability. 
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FIG. 3: Coarsening dynamics in the x-direction (a) and in the 
y-direction (b). Prom bottom to top line-cuts of the height u 
at t = 1000, 2000, 3000, 5350. 



IV. 



ANISOTROPIC COARSENING OF THE 
VICINAL SURFACE 



We integrate numerically (40 1 using a pseudospectral 



method for space discretization and the exponential time 
differencing Runge-Kutta method for time steppin^^. 
Periodic boundary conditions are imposed. The typical 
space grid is 1024^ with Ax = HqL^ — lloLy, and the 
time step At = 0.01 {Iq/^IDCo)T. Resolution and size 
of the simulation are enough to obtain reliable statistical 
data. Throughout we use units = IqL^, ly = loLy, 
to — (Iq/UDCq) T, and uq — Iq U, for lengths, times and 
heights, respectively. 

We show the evolution of the vicinal surface at differ- 
ent times in Fig. [5] The horizontal axis is in the step- 
wise direction y, it gives the shape of the meanders, and 
the vertical one is in the terrace direction x, it follows 
downhill the mean slope, the third axis represent the 
fluctuations of the surface shape u{x,y,t). We observe 
that the height scale, as well as the size of the structures 
steadily increase in time, which are characteristics of a 
coarsening dynamics. The surface develops in time an 
anisotropic pattern with parabolic meanders spanning in 
the y-direction together with complex fluctuations in the 
x-direction. In order to visualize the shape change of 
the surface in both directions we represent in Fig. |3] a 



6X) 

o 



-2 
-4 



jS= 1.06 



log(0 

FIG. 4: Roughness as a function of time. Logarithmic plot 
showing the power law in ui ~ t. 



cut of the height (for example at the edge of the box) at 
fixed y and x. Coarsening is observed along both direc- 
tions, although the coarsening dynamics differs between 
step flow (x-direction) and meanders (y-direction). In 
the y-direction the meanders are composed by a series 
of parabola-like segments reminiscent to the evolution of 
the one dimensional meandering 

instabilitjJ^i] (Fig. M)). 
In particular, at late times, the system instead of tend- 
ing towards a quasi-one-dimensional in-phase state with 
large meanders as would be dictated by the sole linear 
instability, a full two-dimensional state with persisting 
fluctuations remains (c.f. Fig. [2|d)). 

Using data from the time evolution one can compute 
the roughness 



w{t) = ^{u{x,y,tY)-{uy- 

(in our case the mean value {u) = vanishes) it is repre- 
sented in the graph of Fig. [4] in logarithmic scales. We 
can fit the long time roughness by a power law w ^ 
with exponent (3 = 1.06. 

In order to characterize the inherent anisotropic evolu- 
tion of the surface morphology it is convenient to define 
characteristic lengths along both directions, perpendicu- 
lar At, 



J dkdq fc^ 



\UkQ 



and parallel Aj^ 



/ dkdq \ukq\ 
to the steps, 

/ dkdqq^ \ukqf' 



J dkdq \ukq\ 



(41) 



(42) 



where Ukq{t) is the spatial Fourier transform of u{x, y, t). 
Figure [5] shows the behavior of \x,y{t) and the power law 
fits. At long times we obtain exponents of Ux = 0.25 and 
ay — 0.52, for A^; and Xy, respectively. This behavior can 
be correlated with the evolution of the power spectrum 
of the height fluctuations 



P.(fc) 



Idq\ukq{t)\^ 
J dkdq |wfeg(i)|" 



8 




- 



1000 2000 3000 4000 5000 

t (a) 





(b) 



(c) 



FIG. 5: Coarsening lengths as a function of time, (a) The 
characteristic length increases faster in the meandering direc- 



tion, A„ 



1/2 



in red, and fit in (b), than in the step flow 

,1/4 



direction. Ax ~ t in blue, and flt in (c). 



and 

^ jdk\uk,{t)? 
^^^'^^ S dkdq\uu,{t)\^ ' 

as shown in Fig.|6] The spectrum of the terrace- wise fluc- 
tuations is rich at small wavevectors and shift towards 
the long wavelength direction in time; the stepwise fluc- 
tuations are peaked around a well defined wavelength 
(reminiscent of the initial most unstable mode), and as 
the coarsening develops the peak shift towards the long 
wavelength direction. 

It is straightforward to explain the power laws observed 
in the numerical simulation of (40 1; a simple similarity 
argument suffices. Equation (401 satisfies the conserva- 
tion relations, 



— / dxdy u = , 
dt J 

the mean value of u is constant, and 



(43) 



d 
di 



J dxdy ^ J dxdy [{dyu}'^— 



C{d^yuf] (44) 



the amplitude grows by the instability term preferentially 
in regions of str ong gradients; therefore, as in the one- 
dimensional caseSSl the amplitude increase is controlled 
by the balance of the instability with the nonlinear term. 
On the other hand, we note that at long times the dis- 
persive wave terms are dominated by the instability, re- 
laxation and nonlinear terms. Indeed, linear terms give 




FIG. 6: One-dimensional power spectrum (a) in the x vicinal 
slope direction, and (b) in the y step-wise direction. Three 
times are shown: 1000 (blue) 3000 (red) and 5350 (black). 



an increase in the characteristic wavelength x,y ^ t^^^, 
faster than the relaxation length evolution x, y ~ t-^^'^. 
Therefore, the effective long time equation contains es- 
sentially the driving instability term which competes with 
the nonlinear one, and the slowest relaxation term in the 
forth derivative of x: 



dtU ^ -dxXXxU - dyy [u + {OyUY 



Putting the similarity ansatz. 



u{x,y,t)^tf^u[-^^ 



V 



(45) 



(46) 



into ( [45} one obtains the power law exponents P — 1, and 
ax — 1/4, ay = 1/2, close to the numerical results (c.f 
Fig. [5}. We confirm that this scaling is compatible with 
( [44| , where the time derivative and the first two terms of 
the right hand side are of the same order, they increase 
linearly in time, while the two last terms rapidly become 
negligible {{dyyu)"^ ^ const, and {dxyuY ^ t^/^). 

We note that the long time effective dynamics de- 
scribed by Eq. ( [45] ) is independent of the physical param- 
eters. It accepts, as the original equation pO| , particu- 
lar one-dimensional solutions u — u{x,t) or u — u{y,t). 
However, the equation is not separable (the product func- 
tion U{x,y) = Ux{x)Uy{y) is not a solution, because of 
the nonlinear term), showing that the long time behavior 
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.^ittl 



1 2 3 4 5 

FIG. 7: Atomic force microscopy (AFM) image of Si(OOl) 
from experiments of Pascale et al!^ after deposition of 300 
monolayers. The system's size is 5 x Inset shows the 

Fourier power spectrum. 



is essentially two-dimensional. The fact that the charac- 
teristic exponents of the meander coarsening are similar 
to the ones of the in-phase case, explains the reminiscence 
of the observed behavior to the one-dimensional case. In 
addition, as can be inferred from the partial spectra of 
Fig. |6] the time dependence of the characteristic length 
scale of the step flow A^; is related, not to structures of 
Aa; size, but to the steady shift of the spectrum towards 
small wavenumbers. This behavior contrasts with the y- 
direction spectrum, dominated by the characteristic size 
Ay, of the parabolic meanders. Therefore, an appropri- 
ate description of the meander dynamics can be given by 
the coalescence of neighboring parabolic meanders, but 
at variance to the one-dimensional in-phase case, per- 
turbed by the persistent step flow fluctuations. (It is 
worth mentioning, that adding a white noise term to ^ 
does not change its scaling behavior, as we found from 
direct numerical computations.) 



and characteristic lengths parallel or perpendicular to 
the steps. As an interesting consequence an initial out 
of phase perturbation of the straight steps does not tend 
at long times towards the (most unstable) in phase pat- 
tern even if this symmetric mode is the most unstable 
one. Dephasing persists due to the slower coarsening in 
the step flow direction, A^, ~ t^^'^, than the coarsening of 
meanders along the steps direction, Aj, ~ t^/^. 

Although the present model can describe the weak non- 
linear regime of the meandering instability under the con- 
dition that the step flow is stable, it would be important 
to generalize Eq. (40 1 in order to include the step bunch- 



ing instability. The two instabilities can be present simul- 
taneously even in ho moepi taxy, in the case of anisotropic 
diffusion as in Si(001)'22H2] Experimental evidence of the 
persistence of two-dimensional patterns is presented in 
the atomic force image of Fig. [7] where we see a typical 
Si(OOl) vicinal surface grown by molecular beam epitaxy, 
after the deposition of 300 monolayers.!^ jjj particular, 
at this late stage of the surface evolution, the form of 
the Fourier spectrum (inset) is somewhat similar to the 
one of Fig. l6l but with the roles of and ky inverted. 



A possible form of the bunching-meandering instability 
equation is. 



dth = L~ V^lV/ip - AV • {d.^hVh) 



(47) 



where V = {dx, dy) and the linear operator L is in Fourier 
space, of the form 

L = ak^-k'^ + bcf -q^ + ik{Ak^ + Bq^)~ Ck^q^ , (48) 

with a, b constants that control the bunching and 
meandering instability growth rates, respectively, and 
X,A,B.C — const, depending on the physical parame- 
ters, in particular A must be proportional to the flux. 
Indeed, the nonlinear term is a generalization of the 
Fdx{l/hx) term, proportional to the flux, appearing in 
the unstable step flow. Derivation of (47 1 from a mi- 
croscopic model deserves further investigations. More 
generally the method introduced in this paper, valid in 
the weak amplitude, long wavelength approximation, ap- 
plied here to the meandering instability, can be adapted 
to more general cases, notably to the case of heteroepi- 
taxial growth of thin films on vicinal surfaces. 



V. CONCLUSION 

In this paper we devised a method to obtain the evo- 
lution equation of a vicinal surface from the continuum 
limit of the adatoms diffusion microscopic model. We 
demonstrated that the anisotropy inherent to a vicinal 
surface induces different time scales for the amplitude 
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